How to derive and parameterize effective potentials in colloid-polymer 

mixtures 
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Polymer chains in colloid-polymer mixtures can be coarse-grained by replacing them with single soft particles 
interacting via effective polymer-polymer and polymer-colloid pair potentials. Here we describe in detail how 
Ornstein-Zernike inversion techniques, originally developed for atomic and molecular fluids, can be generalized 
to complex fluids and used to derive effective potentials from computer simulations on a microscopic level. In 
particular, we consider polymer solutions for which we derive effective potentials between the centers of mass, and 
also between mid-points or end-points from simulations of self-avoiding walk polymers. In addition, we derive 
effective potentials for polymers near a hard wall or a hard sphere. We emphasize the importance of including both 
structural and thermodynamic information (through sum-rules) from the underlying simulations. In addition we 
develop a simple numerical scheme to optimize the parameterization of the density dependent polymer-polymer, 
polymer-wall and polymer-sphere potentials for dilute and semi-dilute polymer densities, thus opening up the 
possibility of performing large-scale simulations of colloid-polymer mixtures. The methods developed here should 
be applicable to a much wider range effective potentials in complex fluids. 
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I Introduction 

Binary mixtures of colloidal particles and 
non-adsorbing polymers have received renewed 
and growing attention recently, in part because 
they exhibit complex and interesting structure, 



phase 
ologyl 



ill 



facial properties, and rhe- 
Hi30, and in part because they 
are excellent model systems for the study of large 
length and time-scale separations in complex flu- 
ids. Problems with bridging length-scales are im- 
mediately apparent in even the simplest models 
of colloid-polymer mixtures: while the mesoscopic 
colloidal particles can be modeled as hard con- 
vex bodies, the polymers are often treated at the 
microscopic (Kuhn) segment level. Thus, even 
though the average size of the polymer coils may be 
of the same order of magnitude as that of the col- 
loids, the number of degrees of freedom needed to 
model the former may be several orders of magni- 
tude larger than what is needed for the latter. This 
naturally provokes the question: Can the polymers 
also be modeled as single composite particles? In 
fact, this is exactly what was done by Asakura 
and Oosawa (AO) who, iB|-t.hcir classic work on 
colloid-polymer mixturesE3t3, modeled the poly- 
mers as ideal particles with respect to each other, 
and as hard-spheres with respect to the colloids. 
This model is strictly speaking only valid for non- 



interacting polymers, or for interacting polymers 
in the dilute limit, while many interesting phe- 
nomena such as polymer induced phase separation 
take place at finite concentrations of interacting 
polymers. Our ultimate goal, therefore, is to go 
well beyond the AO model and describe non-ideal 
polymers in a good solvent up to semi-dilute den- 
sities. We recently extended the AO concept to 
take into account polymer-polymer interactions, 
first by rather naivelyr&dding a Gaussian repulsion 
between the polymeraij, and then by carrying out 
a much more sophisticated programme which re- 
sulted in density dependent polymer-polymer and 
polymer-wall pair interactionsE£IEj, or, in a com- 
plementary approach, density independent many- 
body interactionstH. 

In this paper we revisit the route to the den- 
sity dependent pair potentials in,, much more de- 
tail than in our previous worUi3tL3. These den- 
sity dependent interactions are derived by invert- 
ing structural information, namely the center of 
mass (CM) radial distribution function g(r) for 
the polymer-polymer interactions, and the wall- 
polymer or sphere-polymer CM density profile p(r) 
for the polymer-wall or polymer-sphere interac- 
tions respectively. We focus here on the consid- 
erable subtleties inherent in these inversions, in 
particular the importance of using thermodynamic 
information or related sum-rules to achieve accu- 
rate potentials. Many of the lessons learned should 
be valid for a wider set of effective potentials. As 
a first example, we derive effective potentials for 
polymer solutions based on mid-point and end- 
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point representations. 

We also revisit the problem of deriving effective 
potentials for a wall-polymer interaction, again 
emphasizing the subtleties involved in the inver- 
sion process. In addition, we show how to use a 
similar method to derive accurate sphere-polymer 
potentials, a key step towards large-scale simula- 
tions of colloid-polymer mixtures. Whereas the 
integral over the potentials holds the key to ac- 
curate thermodynamics in the homogeneous poly- 
mer case, here the relative adsorption, defined as 
the integral over the density profile near a wall or 
sphere, is the key to achieving the correct thermo- 
dynamics. 

We also derive a parameterization scheme for 
the potentials used in the homogeneous and the 
inhomogeneous systems. To achieve this, we use a 
novel Monte-Carlo scheme which should be easily 
adaptable to a wider class of effective potentials. 

The paper is organized as follows: In Sec. II we 
repeat the most important aspects of the models 
and simulation methods we used. In Sec. Ill we 
describe how to derive polymer-polymer, polymer- 
wall, and polymer-sphere effective potentials from 
self-avoiding walk (SAW) polymer simulations, 
and discuss the subtleties involved in the inversion 
procedures. We also derive other representations 
of the polymer-polymer potential, such as the mid- 
point or end-point representation. Section [V con- 



tains the fitting procedures we use to parameterize 
the density dependent potentials. We end with the 
customary concluding remarks. 



II Simulation Models and Methods 

Simple lattice models, such as the self avoiding 
walk (SAW), are widely used to describe flexible, 
uncharged polymers in a good solvent. Their popu- 
larity stems in part from the fact that, even though 
they ignore all chemical details except excluded 
volume and polymer connectivity, these models 
can still reproduce the scaling behavior and many 
other physical properties of athermal polymer solu- 
tionsEilH. In addition, because of their simplicity, 
these models are ideally suited for highly efficient 
numerical simulation algorithms. 

The SAW model consists of a cubic lattice of size 
M on which N chains of L monomers are placed. 
Repulsion between (monomer) segments is built in 
by allowing only one segment per lattice site. The 
monomer concentration is given by c = NL/M, 
while the polymer chain concentration is given by 
p = N/M. The size of the polymers is determined 
from the radius of gyration R g , which, for an iso- 
lated chain, scales as R g ~ L 59 . 

Solutions of flexible polymers in a good solvent 



can be classified as a function of polymer con- 
centration into the dilute, semi-dilute and melt 
regimes, each with different scaling behavior. The 
density or concentration at which, roughly speak- 
ing, the dilute regime crosses over to the semi- 
dilute regime is called the overlap concentration 
or density p* , defined as ^7rp*R g = 1, and cor- 
responding to 1 polymer per unit volume ^%R g . 
In both these regimes, the monomer concentration 
c <C 1, and the only relevant length-scale is R g . 
When, upon increasing the polymer density, the 
monomer concentration becomes appreciable, the 
system crosses over to the so-called melt regime 
where the monomer size becomes an additional 
relevant length-scale. Thus, when modeling the 
semi-dilute regime, it is important to take poly- 
mer chains that are long enough, to ensure that c 
is still small. The monomer density c* at the over- 
lap concentration p* scales roughly as c* w 4L~ 0,8 
for SAW lattice chains (see Ref.t£l), which implies 
that for polymers of finite length, there is only 
a limited semi-dilute regime. We studied poly- 
mers with lengths L = 500 and L — 2000. For 
an isolated L = 500 coil the radius of gyration is 
R g = 16.5 ± 0.02 so that c = 0.027 at p = p* , 
and c = 0.23 at p/p* = 8.70, the highest density 
we study for this length. For L — 2000 chains, 
R g = 37.45 ± 0.04, so that c = 0.009 at p = p* , 
and c = 0.064 at p/p* = 7.04, the highest density 
studied for L — 2000 polymers. This suggests that 
while we do not expect the L — 2000 polymers to 
exhibit any significant corrections to scaling behav- 
ior for the densities studied, the L — 500 polymers 
may begin to deviate slightly from the true semi- 
dilute regime at the highest densities. We estab- 
lished previously that the second virial coefficient 
for two L = 500 polymers is close to the scaling 
lirnitfia. This indicates that for low densities, the 
effective potentials, obtained by the coarse grain- 
ing procedure, are very near the scaling limit where 
properties only depend on R g and not on L. 

The SAW simulations were performed JsjL-em- 
ploying the Monte Carlo pivot algorithmEj-O to- 
gether with translational moves. This simple al- 
gorithm is sufficient to sample the configuration 
space at low concentration. In the semi-dilute 
regime, i.e. for concentrations p > p* , we . use C on- 
figurational Bias Monte Carlo algorithmsEjCJ. 



Ill Deriving effective pair potentials 

The central theme of the coarse graining wn 00 - 
dure advocated here and in previous workE£lll3, is 
to replace each polymer by a single particle, in- 
teracting with the other polymer particles via an 
effective (pair) potential. For complex particles 
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like polymers, there is some freedom in the choice 
of coordinate for the single "particle". We have 
mainly used the CM, but one could also use the 
average monomer, the cndpoints, or the midpoint 
as a monomer based— representation (see e.g. the 
appendix of Refc3 ortZl) . We therefore discuss the 
midpoint and endpoint representations in section 
[II B; a more conmlete analysis is planned for a 



future publicatior 



A CM-CM potentials for homogeneous poly- 
mer solutions 

In the p — > 0, or low density limit, the effec- 
tive potentials can be derived from the logarithm 
of probability of overlap of the CM of two poly- 
mer coilstJ'EZl. Calculating this overlap probability 
involves integrating over the polymer degrees of 
freedom by the Monte Carlo procedures described 
in Sec. II. We sample the configurations of two 
polymers an infinite distance apart with the pivot 
algorithm, and after every 1000 pivot moves we 
accumulate an overlap probability histogram by 
testing for segment overlap as a function of the 
CM distance. The effective potential determined 
in this manner is approximately Gaussian in shape, 
with a range of the order R g for all lengths. When 
the CM's completely overlap ,the potential has 
a maximum of v(r = 0) = 1.88 ± 0.01 for the 
L = 500 polymers and v(r = 0) = 1.82 ± 0.02 
for the L = 2000 polymers, which is very close 
the scaling limit estimate v(r = 0) = 1.80 ± 0.051 
(Since all interactions are of entropic origin, we 
set f3 — l/kgT = 1 throughout this paper). In 
the scaling limit the potentials depend only on R g , 
so that the probability for complete overlap of the 
CM's of two polymers will be independent of their 
length L. 

At finite densities there is no longer a simple 
logarithmic relationship between overlap probabil- 
ities and the effective potenJiaM£lEZl. Instead, we 
use the one-to-one mappingejo between g(r) and 
v(r) to find, for each density, the unique effective 
potential that exactly reproduces the two-body 
correlations. To generate the pair-correlations, we 
performed MC simulations of L = 500 polymers on 
a 240 x 240 x 240 cubic lattice, for several different 
number of polymers ranging from N = 50 (p/p* = 
0.07) to N = 6400 (p/p* = 8.70). For L = 2000 
polymers we used a lattice of 500 x 500 x 500 units, 
and the number of polymers ranged from N = 200 
(p/p* = 0.35) to N = 4000 (p/p* = 7.04). During 
the simulations we collected the radial distribution 
functions between the CM of the polymers. We 
typically needed on the order of 10 7 Monte-Carlo 
moves to achieve sufficient accuracy. 



The effective potentials v(r; p) were constructed 
using the one-component OZ equation—supple- 
mented by the hypernetted chain closureo, which 
is nearliz-. exact for the type of potentials we 
generated. The resulting effective pair interac- 
tions are plotted in Figure | and Fi gure 2 for 
the L — 500 and L = 2000 polymers respectively, 
and show a clear though surprisingly small density 
dependence. This density dependence can be un- 
derstood from the effect .of the density- independent 
many-body interactionstll which, for p/p* < 1, are 
dominated by the three-body interactions. 

The basic principles of the_inversion procedure 
were already described in RefEEl, but here we want 
to point out some important subtleties that must 
be kept in mind when performing such inversions. 

(1) First of all, it is important to generate accu- 
rate g(r)'s. This is illustrated in Figure ||, where 
two different potentials from Figure [I], namely the 
potential at p — and the potential at p/p* = 1.09, 
were used at the same density p/p* = 1.09 to gen- 
erate the radial distribution function g(r). These 
g(r)'s are then compared to the true SAW g(r) 
at that density. Even though the potentials seem 
quite different (see Figure |lj), the differences in 
the radial distribution functions are small, imply- 
ing that the process which generates the g(r) must 
be significantly more accurate than the difference 
between the radial distribution functions in Figure 
U if one is to obtain accurate effective potentials. 

(2) Secondly, because polymers do not have 
a hard-core, one can achieve much higher num- 
ber densities than is normally found for simple 
fluidsEaEj. This also puts extra demands on the 
accuracy of the correlation functions, since at high 
densities very small errors tend to destabilize the 
OZ inversion procedure, making convergence diffi- 
cult to achieve. 

(3) Thirdly, g(r) is only known up to half the 
simulation box size, rendering the inversion prob- 
lem underdetermined. If the potential is known 
for all r, OZ techniques can be used to extend g(r) 
beyond the box-sizetil, but when one is trying to 
find g(r) from simulations with an unknown v(r), 
some assumption for the large r behaviour of v(r) 
must be made. We assume that v(r,p) — be- 
yond half the box size, but this is not necessarily 
obvious at high density, and one needs large simu- 
lation boxes to achieve proper convergence of the 
effective potentials. We used boxes with sides of 
approximately 14i? g . Although it may seem that 
the potentials have vanished already at shorter dis- 
tances, this is still near the minimum length nec- 
essary to ensure that there are no cut-off effects 
in the potentials, especially for higher densities. 
Note that our new Monte-Carlo fitting procedure, 
described in section IV, partially helps overcome 
this problem. 
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An effective potential that correctly reproduces 
pair-correlations should, in principle, also pre- 
dict the correct thermodynamics through the com- 
pressibility equationca But we stress that cor- 
rectly determining the equation of state (EOS) 
Z = 11/ p (here II is the polymer osmotic pressure) 
or other thermodynamic properties from effective 
potentials can also be quite subtle. We illustrate 
this in Figure ^ where we compare a typical poten- 
tial and its fit to a single Gaussian. Although these 
potentials do not seem that different, and gener- 
ate almost identical radial distribution functions 
g(r), they result in pressures which differ by about 
10%. The origin of this difference is clear from the 
inset, where we plot r 2 v(r;p), which is in fact a 
better measure of the relevant differences between 
the potentials than v(r;p) itself. That r 2 v(r;p) is 
a better indicator for the accuracy of the predicted 
pressure follows from the fact that potentials such 
as those shown in Figure .-El— Figure ^, or Figure 
|] lead to mean-field fluidsEE3, so named because 
their equation of state approximately follows the 
simple mean-field form: 

Z = P/p^l + h(0;p)p, (1) 

for a wide range of densities. Here, '0(0; p) is 
the k = component of the Fourier transform of 
the pair interaction, which for spherical symmetric 
functions involves the integral over r 2 v(r; p). 

Keeping these subtleties in mind, we found very 
good agreement between the EOS generated by ap- 
plying the compressibility equation to the effective 
potentials in Figure gand Fi gure B and the EOS of 
the underlying SAW polymer solutions. The com- 
pressibility was calculated with the accurate HNC 
equation. The L = 2000 results almost exactly fall 
onto the L = 2000 EOS computed directly by the 
method of DickmarO. For p/ p* < 5, the L = 500 
EOS calculated with the effective potentials is very 
close to the directly computed one, but at higher 
densities some small differences develop. We at- 
tribute these to the difficulties in achieving high 
accuracy for the tails of the potentials when the 
inversions are performed at these higher densities. 
We expect the L = 500 EOS to be slightly higher 
than the L — 2000 one at higher p/p* because the 
monomer density c is higher for the shorter poly- 
mers, which induces L dependent corrections to 
the scaling limit at the higher densities. Both the 
L = 500 and the L — 2000 EOS, where c remains 
very small at the densities probed, are slightly 
higher than the EOS deijwed by renormalization 
group (RG) calculationsE3a. They approximately 
follow the the correct 11/ cu^ p 13 scaling expected 
for the semi-dilute regimecil up to the highest den- 
sities. In contrast, if one were to use the v(r; p = 0) 
potential at all densities, the EOS would be under- 



estimated and would follow mean-field fluid behav- 
ior with the incorrect H/p ~ p scaling at large p 
instead. So, even though Figure || shows that 
the p = potential results in pair-correlations 
(feC?*) that are similar to the true <?2(?")'s, the effec- 
tive thermodynamics can differ significantly. Since 
the density dependence of the effective pair po- 
tentials was_shown to arise from the many-body 
interactionsEJ, this immediately suggests that the 
corrections to the simple linear p scaling of the 
EOS are due to the three and higher body inter- 
actions. 



B Other representations of potentials for ho- 
mogeneous polymer solutions 

As mentioned at the beginning of this section, 
one could also use other representations for effec- 
tive pair potentials between polymers. For star- 
polymers, for example, the mid-point is a more 
natural coordinate!!!! The / = 2 arm limit of a 
star-polymer would correspond to a normal linear 
polymer, but in the mid-point representation. 

In Figure || we compare the mid-point, end- 
point, and CM representations of the interaction 
between two isolated polymer coils, all calculated 
in the usual way by taking the logarithm of the 
overlap probability for L = 500 SAW polymers on 
a lattice. We also compare the / = 2 arm limit 
of th|S-|Star-polymer interaction derived by Likos 
et alS3, as well as an improvement especially tai- 
lored for linear polymersE! The mid-point and 
end-point representations have a divergence at full 
overlap that scales as lim r _ >0 ^( r ') °c \n(r/R g ), so 
that they appear to be quite different from the 
CM representation. However, because the diver- 
gence is integrable in 3-dimensions. these poten- 
tials still result in mean-field fluidso. This is illus- 
trated in Figure ^| where we plot Anr 2 v{r). The 
divergence has disappeared, and we see once again 
that it is the tails of the potentials, and not the 
small r behavior, which matters for the thermody- 
namics. While in Figure ^ the differences between 
the 2-arm limit of the star-polymer potential and 
the true mid-point potential are quite small (just 
fractions of ksT), Figure |?], together with eq |l|, 
shows that they will generate dramatically differ- 
ent EOS. But even the specially tailored mid-point 
potential, which appears to almost perfectly follow 
the simulations in Figure ^|, differs visibly when 
plotted as r 2 v(r). In contrast, for the SAW sim- 
ulation, the three different representations, mid- 
point, end-point, and CM, should all lead to the 
same EOS, provided that, as-e*pected, the volume 
terms are similar and smalll 10 H 16 l. 

These observations demonstrate a more general 
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point, also mentioned in the previous section: Just 
because a fit to a potential appears to be very close 
to the original potential, does not necessarily mean 
it will generate the correct thermodynamic behav- 
ior; great care must be taken to ensure that the fit 
conserves the right quantities. (Note that the fit 
m RefB was indeed constrained to give the cor- 
rect 2nd virial coefficient, which will result in the 
correct low density thermodynamics.) 

In principle, one can also derive density depen- 
dent potentials for the mid-point or end-point rep- 
resentations, just as was done for the CM repre- 
sentation. In Figure ||, we show the mid-point po- 
tentials which reproduce the mid-point g(r)'s gen- 
erated from direct SAW polymer simulations. We 
used the same HNC inversion techniques used ear- 
lier for the CM representation, and expect similar 
accuracy. In summary then, while there is some 
flexibility in the choice of coordinates, many of 
the lessons learned for the CM representatipii carry 
right through to the other representations^. 



C Colloid-polymer potentials from wall-OZ 
and sphere-OZ inversions. 

So far, we have only considered effective 
polymer-polymer pair potentials, but a full coarse 
grained description of a colloid-polymer mixture, 
or, more generally, of polymers in confined geome- 
try, requires effective potentials between polymers 
and (colloidal) surfaces as well. We focus here on 
non-adsorbing surfaces, and calculate the effective 
interaction between polymers and a hard wall, and 
between polymers and a hard sphere (HS). The for- 
mer is important for such systems as mixtures of 
polymers and platelets, or mixtures of polymers 
and very large spheres; the latter helps provide a 
model for mixtures of spherical colloids and larger 
polymers. 

Near a hard non-adsorbing surface, entropic ef- 
fects create a polymer depletion layer, even if 
the chains themselves are non-interacting. Al- 
though there is no polymer-polymer interaction 
when modeling such ideal chains, there will still 
be a finite wall-polymer potential 4>{z). This in- 
teraction can be found for ideal polymers near 
a wall by the simple result <f>(z) = — \n(p(z)/p), 
where p(z) is the polymer density at a distance z 
from the surface and p is the uniform bulk den- 
sity. For polymers in an end-point representation, 
the interaction reduces to <fi(z) = — ln(erf(z/2i? g )), 
while for polymers in the mid-point representa- 
tion (f>(z) = -21n(erf(z/\/2i? ff ))El A similar re- 
sult should follow in the CM representation, but 
we have not yet succeeded in obtaining an analytic 
form. 



Asakura and OosawaOEj modelled polymers as 
ideal w.r.t. each other, and with a hard sphere po- 
tential <p(z) of range R g w.r.t. spheres or walls. 
Here we extend their work and calculate the (j)(z) 
that exactly reproduces the density profiles p(z) 
measured by direct simulations of SAW polymers 
near walls or spheres. 

First we consider SAW polymers near a hard 
wall. The effective polymer-wall potentials (f>(z; p) 
can be derived from the p{z) by combining the wall 
Ornstein Zernike equation with the HNC closure, 
leading to 

^z-p) = -Hp(z)/p)+p J dr> (p{z')/p-l)c b (\v-v'\). 

(2) 

Here, Cb(r) is the bulk polymer direct correlation 
function. The details of this inversion procedure 
are given in Ref.ta. In contrast to the homoge- 
neous case, the inversion involves only one itera- 
tion, because the bulk polymer-polymer potentials 
and hence Q,(r) are already known. 

Using the same explicit SAW polymer model 
as in Section II, we performed MC simulations of 
polymers of length L = 500 on a lattice of size 
M = 240 x 240 x 240 with hard planar walls at 
x = and x = 240. The polymer segments were 
not allowed to penetrate the walls. The simula- 
tions were done for N = 0, 100, 200, 400, 800, 1600, 
and 3200. During each simulation, we computed 
the density profiles p{z), where z denotes the dis- 
tance between the polymer CM and the wall. The 
normalized profiles h{z) = p(z)j p — 1, for differ- 
ent bulk concentrations p, are shown in Figure ^. 
The wall-polymer potentials 4>{z\ p) were obtained 
using the wall-OZ-HNC inversion procedure and 
are plotted in Figure [l0| The wall is essentially 
still hard, but is now softened by an additional 
exponentially decaying repulsion. The density de- 
pendence is more pronounced than in case of the 
effective potentials for bulk polymer solutions. 

The potential for a polymer coil interacting with 
a colloidal HS of diameter a can be found from 
the density profile p(r) around the sphere. Here 
r is the distance from the center of the colloidal 
HS. The inversion procedure is very similar to that 
of the wall-polymer potential geometry. For a bi- 
nary mixture of two components labeled c and p, 
in which colloidal component c is infinitely dilute 
(pc — > 0), the binary OZ equations decouple, and 
reduce to 

h cp {r) = c cp (r) + p J dr'h cp (r')c b (\r ~ r'|). (3) 

Where h cp (r) = p(r)/p — 1, and Cf,(r) is again 
the bulk polymer direct correlation function which 
is determined independently. The above equation 
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can be combined with the HNC closure to yield 
<f>(r; P ) = - H P {r)/p) +pj dr'h cp (r')c b (\r - r'|). 

which can be solved in one step by Fourier trans- 
formation. Figure |ll| shows the polymer density 
profiles h cp (r) around a single colloidal HS par- 
ticle with a diameter of a = 2R a . The effective 
potentials are plotted in Figure [L2[ The interac- 
tion appears to be somewhat softer than for the 
planar wall case, at least in the sense that it is still 
finite at distances that correspond to the CM be- 
ing "inside" the hard colloidal particle. This hap- 
pens more readily for smaller size-ratios a/(2R g ) 
because the polymers can deform around the col- 
loids. Note that this penetration into the colloidal 
HS will not occur in the mid- or end-point repre- 
sentations. 



D Connection to scaling theory approaches 

Polymers in the semi-dilute regime have tradi- 
tionally been very successfully understood with 
scaling theories^. There the fundamental unit is 
no longer the radius of gjftajbion R g at zero density, 
but the blob-size £(/9)E3c3, which appears both 
in the EOS, where ILj ~ £(p) -3 , and in the rel- 
ative adsorption T/p ~ £(p)- For the homoge- 
neous polymer solutions, it is not clear from either 
the g(r) or from the effective potentials where the 
blob-size comes in. The reason the effective poten- 
tials correctly reproduce the EOS is because they 
correctly reproduce the pair-correlations, which, 
through the compressibility theoremo, give the 
correct thermodynamics. This blob-size only ap- 
pears in the potentials if one attempts to param- 
eterize them for p/p* > 2 (roughly speaking the 
density at which the scaling sets in), since the in- 
tegral over r 2 v(r) should scale as p~ 2 £{p)~ 3 - Since 
we only parameterize the potentials for lower den- 
sities, this scaling does not apply to the functional 
form we use. 

For polymers near a wall the width of density 
profile of the monomers scales with £(/?) in the 
semi-dilute regimeEl. A similar scaling is expected 
to hold for the width of the CM density profile. 
So even though £(p) is not immediately evident in 
the effective wall-polymer potentials, it is present 
in the induced density profiles. Nevertheless, many 
questions as to the direct relationship between the 
scaling theories and our soft-colloid approach re- 
main to be worked out. Partially for that reason, 
some care must be taken when applying the soft- 
colloid picture deep in the semi-dilute regime. 

The semi-dilute regime scaling theories are not 
normally valid in the dilute regime, where we ex- 



pect our soft-colloid approach to be robust. 

IV Fitting effective pair potentials 

A CM-CM potentials for bulk polymer solu- 
tions 

When applying effective potentials in coarse 
grained simulations or in theoretical analysis, it is 
convenient to have an explicit expression or param- 
eterization of the potentials. In RefJl3 we fitted the 
effective pair potentials we obtained by inversion 
of the g(r) to a single Gaussian and described the 
remainder by means of a spline fit. However, for 
practical use in simulation or for a theoretical anal- 
ysis, one would like to use an analytic expression 
without the complication of a multi node spline 
fit. Moreover, it would be very useful to be able 
to explicitly model the density dependence of the 
potentials as well. Since we are dealing with an 
approximately Gaussian form, the simplest analyt- 
ical expression is a sum of Gaussians with density 
dependent coefficients: 

n 

v f {r-~p) = Y j ai{~p)z- {x,bd ~ P)) \ (5) 

i=0 

where we introduced p = p/p* for clarity. Here, 
the Gaussians are centered at r = 0, because the 
maximum of the potential is located at the ori- 
gin. The density dependence comes in through 
the coefficients a,i(p) and bi(p). In first instance, 
we take the coefficients to be linear in density, i.e. 
a,i(p) = a i0 + CLnp and bi(p) = b i0 + b a p- 

Note that an expression of this form cannot de- 
scribe the slightly negative tails in the potential 
that were found in Ref.EEl. However, as was argued 
in the previous section, the structure is not very 
sensitive to small changes in v(r;p). A more im- 
portant aspect of the fitting procedure is to make 
sure that the fitted potentials lead to the correct 
EOS. To do this we make use of the accurate 
mean field approximation in eq [l| Because this 
is such a good approximation, keeping the k = 
Fourier component v(0; p) = 4?r j r v(r; p) for the 
fitted potential Vf(r;p) equal to that of the origi- 
nal potential v(r; p) ensures that the pressures will 
also be (almost) equal. A fitting procedure should 
therefore include the constraint 

n 

0(0; p) = 0/(0; p) = J £ a t (p) (h(p)f , (6) 

i=0 

which results in one of the coefficients being fixed 
by the constraint. The v(0;p) itself is a smooth 
function of the bulk density p; we fit it to a cubic 
polynomial. 
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The inverted potentials v(r,p) were fitted to 
eq H by applying a standard least squares non- 
linear fitting procedures including the above con- 
straint. To obtain a good fit of all the potentials 
for p/p* < 2 we needed at least 3 Gaussian terms, 
particularly to correctly describe Vf(r = 0; p). For 
higher densities, the linear density dependence in 
eq ^| broke down. This can partly be under- 
stood from the fact that the potential at r = in- 
creases almost linearly for p/p* < 2 but decreases 
for higher concentration (see inset of Figure |l|). 
Clearly, a linear fit in the density cannot cope with 
this behavior. Higher order polynomials for the co- 
efficients were not very successful either. A better 
option would be to fit the potentials for high den- 
sities independently, and then ensure continuity at 
the crossover. 

After using the non-linear fitting procedure to 
directly fit the potentials, the best fit coefficients 
in eq @ can be further fine-tuned by minimizing 
the difference between the radial distribution func- 
tion g(r) of the SAW polymers and the gf(r) gen- 
erated from the fitted potential vj(r). This min- 
imization procedure works as follows. For a cer- 
tain density p — p/p*, we first calculated gj(r) for 
Vf(r; p) with the best fit coefficients using the HNC 
approximation. We then compare this <7/(r) to the 
original radial distribution function gsA\v{ r )i mea- 
sured in the explicit SAW simulation at the same 
density, by determining the least square difference: 
A (p) = J(9f( r ) ~ 9SAw{r)) 2 r 2 dr. We performed 
this calculation for every density pi we have SAW 
data for, summing up to A to tai = J2i Sub- 
sequently, we changed one of the coefficients by a 
small amount and determined the difference Atotai 
again. If the difference between the fitted g/(r) 
and the SAW data gsAw becomes smaller, the new 
coefficient was accepted, otherwise the old coef- 
ficient was restored. These steps were repeated, 
until the difference reached a minimum value and 
changes in the coefficients were no longer accepted. 
The final coefficients only represent a local mini- 
mum, so this Monte Carlo type minimization pro- 
cedure cannot be applied without a good initial 
estimate for the fit coefficients.. On the contrary, 
it depends strongly on the quality of the original 
non-linear fit described above. Note that by using 
this optimization we skip one step, namely the in- 
version of the original g{r). If a small error was 
made in the inversion, for example by our assump- 
tion that v(r) =0 for r > r cuto ff, it should be 
corrected during the optimization. However, one 
must keep in mind that the minimization proce- 
dure only finds a local minimum. Larger errors 
may not be corrected. 

We applied the above fitting procedure to the 
bulk polymer results for L — 500 and L = 2000 . 
The coefficients are given in Table | and Table O, 



respectively. One of the coefficients, 63 (p), is fixed 
by the constraint of eq H: 



i=l 

(7) 

Note that this implies that in contrast to the other 
coefficients, 63 (/5) is no longer linear in density. 
The best fit coefficients for both lengths are not 
far from each other, again indicating that the re- 
sults are close to the scaling limit. 

To illustrate the importance of explicitly includ- 
ing the constraint of eq ^, we applied the same 
non-linear fitting procedure, including the Monte- 
Carlo optimization, but without the constraint. 
In Figure (|l3|) we compare v(0;p), which is a 
good measure for the EOS, for the constrained 
and the unconstrained case. Clearly, including the 
constraint is essential for reproducing the correct 
thermodynamic behavior. At this point we would 
like tp-.return to a statement made in an earlier 
papenlj, where we emphasized the importance of 
the very small negative tails we found from our in- 
version procedure. Because the tails are so small, 
their effect on the radial distribution functions is 
almost imperceptible. It is only when they are 
multiplied by r 2 that they become important for 
the EOS. Our current fitting procedure ignores the 
negative tails, but because we force the sum-rule 
on r 2 v(r), the thermodynamics are still correctly 
reproduced. Therefore our previous statements on 
the importance of the tails should be tempered. It 
is actually the sum-rule on the pair-potential which 
is critical to achieving the correct thermodynam- 
ics. 



B Colloid-polymer potentials 

A similar fitting procedure as in the previous 
section can be applied to the effective wall-polymer 
potentials. In Ref.E3 we showed that the effective 
wall-polymer potential (f>(z, p) can be fitted reason- 
ably by the cubic exponential function 

4>f(z, p) = a (p) exp [ai(j5)z + a 2 {p)z 2 + a 3 (p)z 3 .] 

(8) 

As before, the density dependence can be intro- 
duced through the coefficients a,i(p) = + anp + 
aiip 1 . Here we chose a quadratic density depen- 
dence, because a linear expression yielded a poor 
fit. We fitted the potentials <fi(z; p) using the non 
linear fit method of the previous section but with- 
out any constraint. Subsequently, the coefficients 
were optimized using the minimization procedure 
mentioned above. The optimized best fit coeffi- 
cients are given in Table [II, and the hf{z) gener- 
ated by these fitted potentials are shown in Figure 
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|9j. The fit seems to be reliable for p < p* , but, 
unfortunately, begins to break down in the semi- 
dilute regime. For the highest densities, the repro- 
duced hf(z) show more structure than the den- 
sity profile obtained from the SAW simulations. 
Clearly, eq || cannot completely capture all the 
wall-polymer effects in the semi-dilute regime. The 
functional form of eq ^ is rather ad hoc, and the 
fit would probably be better if one had access to a 
more accurate analytical expression based on phys- 
ical arguments. 

Similarly to the bulk polymer case, one needs 
to compare surface thermodynamic properties to 
ensure the quality of the fit. We focus on the poly- 
mer adsorption at the wall, T, defined as the par- 
tial derivative of the excess grand potential Q, ex per 
unit area with respect to the chemical potential /i, 



_ d(Q ex /A) _ 
dp P 



h(z)dz, 



(9) 



where A is the surface area. In Figure |lj we show 
that the adsorption T is indeed well represented 
by the fits. So, even though for high densities the 
structure next to the wall is not well described by 
eq ^, the integral over h(z) is still accurately rep- 
resented. Since the adsorption T is given by the 
integral over h(z), instead of an integral weighted 
by z 2 , the tails of the potentials are less impor- 
tant to the correct thermodynamics than in the 
case of the bulk polymer potentials. In fact, if the 
adsorption T and the EOS are known as a func- 
tion of the density, then the surface tension of the 
polymer solution cae-be calculated with the Gibbs 
adsorption equationEa. Therefore, since the fitted 
wall-polymer potentials give the correct T, and the 
polymer-polymer potentials give the correct TLb/p, 
when used together they should provide an accu- 
rate representation of surface tensions and related 
surface thermodynamic quantities. 

In the case of a spherical particle, the fit to the 
sphere-polymer potential is more straightforward. 
Once again we can use a sum of Gaussians to esti- 
mate the potential 



4>f{r) = J2 a i(p) e ~ 



((r- Cl (p))/6,(p)) 2 



(10) 



i=0 



Here, because of the size difference between the col- 
loidal particle and the polymer coil, we allow the 
Gaussians to be off center. In the non-linear fit for 
the a = 2R g data, we set n — 2 and c\{p) = 0, and 
assumed a linear density dependence of the coef- 
ficients. Subsequently, the coefficients were opti- 
mized using the minimization procedure. The re- 
sults are given in Figure [ll] and Table |v|. The 
hf(r) reproduces the measured h(r) quite well. To 
check the quality of the fit we compared in Figure 



|l5| the adsorption of the SAW polymers around 
the sphere with the optimized potentials. Here, 
in contrast to the planar wall case, we do have to 
integrate over h(r)r 2 . We would therefore expect 
that the deviations are larger than for the wall. 
However, although the agreement is not as good 
as in the wall case, the adsorptions are still well 
represented, thus giving us confidence that the fits 
are adequate. 



V Conclusion 

In this paper we have shown in detail how to use 
OZ techniques to derive effective potentials poten- 
tials for polymer solutions from direct simulations 
of SAW polymers. Just as many subtleties were 
found in the original application of OZ inversion 
techniques to simple atomic and molecular fluids, 
so we find that great care must be taken to cor- 
rectly invert and fit our effective potentials. 

We found that g(r) is not very sensitive to differ- 
ences in v(r), which is very similar to the situation 
for atomic and molecular fluidscj. This insensitiv- 
ity places strong demands on the accuracy of the 
original g(r)'s needed as input for our inversion 
procedures. 

For the polymer-polymer potentials, it is cru- 
cial that the integral over r 2 v(r; p) is correctly rep- 
resented. Seemingly very small differences in the 
tails of v(r), which in turn result in almost imper- 
ceptible changes in g(r), can nevertheless result in 
large differences in the EOS. This principle should 
hold not only for linear polymers in the CM, end- 
point or mid-point representations, but alsqJpr low 
arm-number star polymersE£l, dendrimersEa, and 
other mean- field fluids. 

We also derived density-dependent wall-polymer 
and sphere-polymer potentials by directly invert- 
ing the one-particle density profile p(r) calculated 
by direct simulations of L = 500 SAW polymers. 
Here, the important thermodynamic parameter is 
the adsorption T, which is quite well reproduced 
by our fitted potentials. Since our effective po- 
tentials provide a good representation of the EOS 
and of the adsorption, they should lead to an ac- 
curate representation of the surface tension and 
other related surface thermodynamic properties of 
polymer-colloid mixtures. 

And finally, while one might think that most of 
the hard work is done once effective potentials have 
been inverted from direct simulations, fitting these 
potentials for their use in large-scale simulations 
of colloid-polymer mixtures is far from trivial. We 
showed how to use a numerical optimization pro- 
cedure to ensure the accuracy of our fits. Because 
this optimization procedure skips the direct inver- 
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sion step, it can remove residual errors in the in- 
versions, guarantee that the fits conserve the right 
sum-rules, and lead to the correct thermodynam- 
ics. We emphasize that these conclusions should 
hold for a much broader class of effective poten- 
tials than the ones we discussed here. ___ 
The ultimate goal of our research projectrJiljO 
is to model large-scale mixtures of polymers and 
colloids. The coarse-graining of the polymers from 
the "microscopic" SAW chains to single composite 
entities interacting via effective density dependent 
pair potentials is a crucial step toward that goal. 
The next step will be to use the accurate fits de- 
rived in this paper to perform direct simulations 
of many spherical colloids interacting with many 
polymers. 
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Note added in proof We discovered a small mis- 
take in our calculations of the radius of gyration. 
The value for radius of gyration mentioned in 
the text for L = 500 is in reality R g = 16.83 in- 



40 



stead of R„ 



16.5. and R„ = 38.4 for L=2000 



instead of R g = 37.45. This means that the dis- 
tances in terms of R g have to be rescaled by a 
factor of 0.9804 in case of L = 500 and by a 



factor of 0.977 in case of L = 2000. All the re- 
duced densities have to be rescaled by a factor 
of 1.0612 for L = 500 and by 1.069 for L = 2000. 
These small adjustments are only relevant when 
comparing to the theoretical results of the RG. 
However, they do not change the validity of the 
approach nor any of the conclusions. We thank 
Vincent Krackoviack for pointing this out to us. 
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TABLE I: Best fit coefficients for v(r;p) as defined 
in eq 5 for L = 500 polymers. Coefficients for bz(p) 
are fixed by the constraint eq 7. The v(0; p) in eq 
7 is approximated by cubic polynomial v(0;p)/An — 
1.2902 + 0.28132p + 0.13676L5 2 - 0.040892p 3 
i dio an bio foil 

1 1.47409 -0.07689 0.981368 -0.056808 

2 -0.232096 0.031321 0.42123 -0.026278 

3 0.638974 0.24193 - 



TABLE II: Best fit coefficients for v(r;p) as defined 
in eq 5 for L — 2000 polymers. Coefficients for bz(p) 
are fixed by the constraint eq 7. The v(0; p) in eq 
7 is approximated by cubic polynomial v{Q;p)/A-n = 
1.245 + 0.3564275p - 0.02443297/5 2 + 0.0018028/5 3 

i diO dil foiO foil 

1 1.41808 -0.081969 0.979493 -0.057796 

2 -0.224377 0.030647 0.440907 -0.024327 

3 0.630219 0.211378 - 



TABLE III: Best fit coefficients for <j>(z\ p) as defined 
in eq 8 for L — 500 polymers next to a wall. 

i an-, an aa 

62.7242 56.4595 -29.92825 

1 -6.40938 -3.88795 2.044202 

2 2.50812 5.156190 -2.13356 

3 -0.69042 -1.55191 0.59725 
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FIG. 2: The effective polymer CM pair potential 
v(r; p) for L — 2000. Inset: A comparison with Figure 
^ for p/p* — (solid lines), shows that we are very 
near the scaling limit. This is similar for finite density; 
e.g. the effective potential for p/p* = 1.76, L — 2000 
(dotted line) is very close to the p/p* — 2.18, L = 500 
potential (dashed line). 



TABLE IV: Best fit coefficients for <f>(r; p) as defined 
in eq 10 for L = 500 polymers around a sphere of 
diameter a — 2R g . 

i 1 



a i0 5.5610 1.8477 

an -0.8042 1.4759 

fo i0 0.7751 1.2720 

foil -0.1151 0.1052 

Cio 0.4082 0.0 

Cii 0.14104 0.0 
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FIG. 3: Comparison of g(r)'s generated at density 
p = p* by the low-density potential v(r; p — 0) and the 
higher density potential v(r;p/p* — 2.18), compared 
to the true g(r) at that density. Note how small the 
differences are. 




FIG. 4: Comparison between the effective potential 
v(r; p = p*) (solid line) and its fit to a single Gaussian 
(dotted line). Although the fit may appear to be quite 
good to the eye, the pressures obtained from these po- 
tentials differ by about 10%. The reason for this is il- 
lustrated in the inset where the potential is multiplied 
by r 2 , and the differences become more visible. 




FIG. 5: Linear plot of the EOS Z = 11/ p for L = 500 
and L — 2000 polymers, here determined from the 
effective potentials through the compressibility route 
and by direct simulations of the underlying SAW poly- 
mer system. The slight deviations for the L = 500 case 
at the higher densities are most likely due to small in- 
accuracies in the inversion procedure used to generate 
the effectivensptentials. We also show the EOS derived 
from an RGEj approach. Using only the lim p ^o form 
of v(r;p/p*) strongly underestimates the EOS. 




FIG. 6: Endpoint, midpoint and CM representations 
of the interaction v(r) between two isolated polymers. 
Also included apeitworfits to the mid-point potentials 
taken from RefsE3 andE3, as explained in the text. 
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FIG. 7: Endpoint, midpoint and CM representations 
of the interaction 4-irr 2 v(r) between two isolated poly- 
mers. Plotting the potentials in this way accentuates 
the differences. Note in particular ifepw poor the / = 2 
limit of the star-polymer potential^ performs. 




FIG. 8: Midpoint v(r) that exactly reproduces the 
midpoint-midpoint g(r) for different bulk polymer den- 
sities. Inset: 4nr 2 v(r) shows that the more significant 
change is in the tails of the potentials. 
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FIG. 9: Measured h(z) (solid lines) next to a wall 
obtained from SAW simulation compared with the 
hf(z) which follows from the optimized fitted po- 
tential (f>f(z;p) (dashed lines). From left to right 
the curves correspond to bulk polymer densities 
p/p* = 2.27,1.16,0.59,0.30,0.16,0.08 and 0, respec- 
tively. Note that the densities differ slightly from the 
corresponding bulk density due to the depletion at the 
wall. 




FIG. 10: The wall-polymer potential 4>(z\ p) as ob- 
tained from the inversion of density profile p{z). 
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FIG. 11: The polymer density profile h(r) = p(r)/p 
around a colloidal HS of diameter a — 2R g . From 
left to right the curves correspond to a bulk polymer 
density p/p* = 2.18, 1.09, 0.54, 0.27, 0.14 and 0, respec- 
tively. The solid lines represent the SAW simulation re- 
sults, whereas the dashed lines correspond to the h(r) 
that results from an optimized fit to the effective po- 
tential (j)(r; p). 




FIG. 12: The colloid-polymer effective potential <j>(r) 
for a colloidal diameter of o2 = R g . 
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FIG. 13: When the explicit constraint of eq ^ is 
not included, the value of v(k = 0; p) — 4ir J r 2 v{r; p) 
(dotted line) begins to deviate significantly from the 
correct value (solid line). This will have an important 
effect on the EOS. 
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FIG. 14: The relative adsorption T/p (in units of 
R^ 2 ) given by direct SAW simulations of L = 500 
polymers next to wall (open circles), compares well to 
the adsorption calculated from the fitted wall-polymer 
potentials (open squares). 
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FIG. 15: The relative adsorption Fp (in units ol -R,^ 2 ) 
given by direct SAW simulations of L = 500 polymers 
around a colloidal sphere of diameter a — 2R g (open 
circles), compared to the adsorption from the fitted 
colloid-polymer potentials (open squares). Note that 
the agreement is not as good as in the polymer wall 
case, because of the integration over h(r)r 2 



